!--------------------------------------------------------
! $Id: x2y2model.shar,v 1.2 2011/05/14 21:34:35 daw Exp $
! Description:
! Author:  Murray S. Daw <daw@clemson.edu>
! Created: 
!--------------------------------------------------------

module rannums
  IMPLICIT NONE

  INTEGER, PARAMETER, private :: K4B = selected_int_kind(9)
  integer, parameter, private :: r8 = selected_real_kind(10,50)


contains
! from numerical recipes (for F90)
  FUNCTION ran(idum)
    IMPLICIT NONE
    INTEGER(K4B), INTENT(INOUT) :: idum
    REAL :: ran
    INTEGER(K4B), PARAMETER :: IA=16807,IM=2147483647,IQ=127773,IR=2836
    REAL, SAVE :: am
    INTEGER(K4B), SAVE :: ix=-1,iy=-1,k
    if (idum <= 0 .or. iy < 0) then
       am=nearest(1.0,-1.0)/IM
       iy=ior(ieor(888889999,abs(idum)),1)
       ix=ieor(777755555,abs(idum))
       idum=abs(idum)+1
    end if
    ix=ieor(ix,ishft(ix,13))
    ix=ieor(ix,ishft(ix,-17))
    ix=ieor(ix,ishft(ix,5))
    k=iy/IQ
    iy=IA*(iy-k*IQ)-IR*k
    if (iy < 0) iy=iy+IM
    ran=am*ior(iand(IM,ieor(ix,iy)),1)
  END FUNCTION ran

! 
  function rgauss(idum)
! P(x) = Exp[-x^2/2]/Sqrt[2 pi]
! <x^2> = 1
! <x^4> = 3
    implicit none
    real(r8) :: rgauss
    INTEGER(K4B), INTENT(INOUT) :: idum
    real(r8), parameter :: pi = 3.1415926535897932385
    real(r8) :: a1,a2

10  continue
    a1 = ran(idum)
    if (a1.eq.0.) goto 10
    a2 = ran(idum)
    rgauss = sqrt(-2.*log(a1))*cos(2.*pi*a2)
    return
  end function rgauss

function normal(mean,sigma,idum) !returns a normal distribution 
        implicit none
        integer, parameter:: b8 = selected_real_kind(14) 
        real(b8), parameter :: pi = 3.141592653589793239_b8 
        integer gene_size   
        INTEGER(K4B), INTENT(INOUT) :: idum
        real(b8) normal,tmp 
        real(b8) mean,sigma 
        integer flag 
        real(b8) fac,gsave,rsq,r1,r2 
        save flag,gsave 
        data flag /0/ 
        if (flag.eq.0) then 
        rsq=2.0_b8 
            do while(rsq.ge.1.0_b8.or.rsq.eq.0.0_b8) ! new from for do 
                r1=2.0_b8*ran(idum)-1.0_b8 
                r2=2.0_b8*ran(idum)-1.0_b8 
                rsq=r1*r1+r2*r2 
            enddo 
            fac=sqrt(-2.0_b8*log(rsq)/rsq) 
            gsave=r1*fac 
            tmp=r2*fac 
            flag=1 
        else 
            tmp=gsave 
            flag=0 
        endif 
        normal=tmp*sigma+mean 
        return 
    end function normal
 
end module rannums

      
